Vamos a hacer EDA (Exploratory Data Analysis) de datos municipios de Brasil de 2019 de tres fuentes:
Vamos a suponer que nuestro objetivo final es predecir el PIB a nivel municipio. Esto podría ser útil, por ejemplo, si sabemos que para algunos municipios este dato es poco confiable, mientras para otros no. También si es un dato que se publica con rezago con respecto a otros que están disponibles más rápido.
WARNING: esto es solamente un escenario hipotético. Si de verdad quisiéramos hacer este ejercicio deberíamos analizar muchos aspectos, como:
WARNING 2: en este análisis nuestro objetivo es explorar y entender los datos, formular hipótesis (por ejemplo, acerca de qué features pueden ser importantes), etc.. Esto quiere decir las transformaciones que hagamos acá no necesariamente vayan a ser las mismas que apliquemos cuando entrenemos un modelo.
WARNING 3: el siguiente es un EDA de entre muchos posibles. ¡Seguramente no sea el mejor porque es para dar clase!
library(tidyverse)
library(janitor)
library(skimr)
library(GGally)
library(ggpubr)
library(hexbin)
library(glue)
library(scales)
# library(kableExtra) # dont load it because it breaks skimr
# install.packages("remotes")
# remotes::install_github("rpkgs/gg.layers")
library(gg.layers)
Levantamos los tres datasets:
file_bolsa = "data/working/bolsafamilia_municipios_2019.csv"
file_ibge = "data/working/ibge_municipios_2019.csv"
file_varios = "data/working/varios_municipios_2019.csv"
df_bolsa = read_csv(file_bolsa, show_col_types=F)
df_ibge = read_csv(file_ibge, show_col_types=F)
df_varios = read_csv(file_varios, show_col_types=F)
Inspeccionamos algunas características básicas de los data.frames:
glimpse(df_bolsa)
## Rows: 5,570
## Columns: 4
## $ uf <chr> "GO", "MG", "GO", "MG", "PA", "CE", "BA", "BA",~
## $ num_beneficiarios_mean <dbl> 995.2500, 273.6667, 685.8333, 539.0000, 34341.5~
## $ suma_valor <dbl> 1973921, 461869, 1164947, 866639, 84271320, 513~
## $ municipio <chr> "abadia_de_goias", "abadia_dos_dourados", "abad~
head(df_bolsa)
## # A tibble: 6 x 4
## uf num_beneficiarios_mean suma_valor municipio
## <chr> <dbl> <dbl> <chr>
## 1 GO 995. 1973921 abadia_de_goias
## 2 MG 274. 461869 abadia_dos_dourados
## 3 GO 686. 1164947 abadiania
## 4 MG 539 866639 abaete
## 5 PA 34342. 84271320 abaetetuba
## 6 CE 1852. 5131327 abaiara
UF (Unidad Federativa) indica el estado.num_beneficiarios_mean es el número de beneficiarios
promedio a lo largo del año.suma_valor es la suma de los pagos totales en el
año.glimpse(df_ibge)
## Rows: 5,570
## Columns: 14
## $ nome_da_grande_regiao <chr> "Centro-oeste", "Sudeste", "Ce~
## $ codigo_uf <dbl> 52, 31, 52, 31, 15, 23, 29, 29~
## $ nome_da_unidade_da_federacao <chr> "Goiás", "Minas Gerais", "Goiá~
## $ codigo_municipio <dbl> 5200050, 3100104, 5200100, 310~
## $ semiarido <chr> "Não", "Não", "Não", "Não", "N~
## $ hierarquia_urbana_principais_categorias <chr> "Metrópole", "Centro Local", "~
## $ vab_agro <dbl> 5536.665, 32756.932, 56244.012~
## $ vab_industria <dbl> 37221.827, 12233.404, 24158.34~
## $ pib <dbl> 202315.97, 127581.83, 346767.5~
## $ pib_percapita <dbl> 23061.21, 18254.66, 17302.04, ~
## $ atividade_maior_vab <chr> "Demais serviços", "Administra~
## $ grado_urbanizacion <dbl> 96.93234, 51.44562, 63.75355, ~
## $ tipo_urbano <chr> "Urbano", "RuralAdjacente", "R~
## $ municipio <chr> "abadia_de_goias", "abadia_dos~
glimpse(df_varios)
## Rows: 5,596
## Columns: 9
## $ area_geografica_2019 <dbl> 147.256, 881.064, 1045.127, 1817~
## $ despesa_ciencia_tecnologia_2019 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## $ despesa_educacao_cultura_2019 <dbl> 13291163, 6061623, 11262623, 124~
## $ despesa_pessoal_encargos_sociais_2019 <dbl> 18676874, 11345153, 25495707, 26~
## $ despesa_transporte_2019 <dbl> 2192450.55, 1294621.67, 946773.6~
## $ exportacoes_2019 <dbl> 538966, NA, NA, NA, 70309528, NA~
## $ importacoes_2019 <dbl> 8894, NA, NA, NA, 42478, NA, NA,~
## $ populacao_2019 <dbl> 8773, 6989, 20042, 23237, 157698~
## $ municipio <chr> "abadia_de_goias", "abadia_dos_d~
despesa_pessoal_encargos_sociais indica gasto en
personal y cargas sociales.municipio ya está normalizado
entre los datasets.Vamos a hacer un join de los datasets usando
municipio como clave (ya sabemos que el campo está
normalizado). Antes de hacerlo es conveniente revisar si hay
observaciones duplicadas o si algún dataset tiene observaciones
(municipios) faltantes.
dups_bolsa = df_bolsa %>% get_dupes("municipio")
## No duplicate combinations found of: municipio
dups_ibge = df_ibge %>% get_dupes("municipio")
## No duplicate combinations found of: municipio
dups_varios = df_varios %>% get_dupes("municipio")
## No duplicate combinations found of: municipio
munis_bolsa = df_bolsa %>% pull(municipio) %>% unique()
munis_ibge = df_ibge %>% pull(municipio) %>% unique()
munis_varios = df_varios %>% pull(municipio) %>% unique()
cat(
length(munis_bolsa),
length(munis_ibge),
length(munis_varios)
)
## 5570 5570 5596
¿Los de df_bolsa y df_ibge son los
mismos?
# intersect(munis_bolsa, munis_ibge)
length(intersect(munis_bolsa, munis_ibge))
## [1] 5570
¿Cuáles de df_varios no están en
df_bolsa/df_ibge?
setdiff(munis_varios, munis_ibge)
## [1] "assunguy_de_cima" "augusto_severo"
## [3] "barcellos" "cimbres"
## [5] "cococi" "conchas_2"
## [7] "entre_rios_3" "mecejana"
## [9] "monsaras" "muribeca_2"
## [11] "nossa_senhora_do_o_de_goyanna" "olivenca_2"
## [13] "palmira" "ponte_do_itabapoana"
## [15] "porto_de_cima" "riacho"
## [17] "santa_cristina_do_pinhal" "santa_izabel"
## [19] "santo_amaro_2" "sao_goncalo_2"
## [21] "sao_joao_marcos" "sao_miguel_2"
## [23] "sao_sebastiao_do_parahyba" "trancoso"
## [25] "vertentes_2" "vila_franca"
## [27] "villa_verde"
¿Cuáles de df_bolsa/df_ibge no están en
df_varios?
setdiff(munis_ibge, munis_varios)
## [1] "campo_grande_3"
¿Cuáles están en todos los dfs?
# aplica la funcion intersect recursivamente sobre la lista
munis_en_comun = Reduce(
intersect, list(munis_bolsa, munis_ibge, munis_varios))
length(munis_en_comun)
## [1] 5569
# Reduce(fn, elementos)
Hacemos el join. Es importante saber qué observaciones estamos conservando o eliminando en cada paso.
df = df_bolsa %>%
full_join(df_ibge, by="municipio") %>% # 5570 municipios
left_join(df_varios, by="municipio")
# excluimos los que estan en varios que no estan en los otros dos (27 municipios)
Renombramos algunas columnas por comodidad:
df = df %>%
rename(
nome_regiao = nome_da_grande_regiao
,nome_uf = nome_da_unidade_da_federacao
,hierarquia = hierarquia_urbana_principais_categorias
,tipologia = tipo_urbano
)
Para normalizar o limpiar los nombres de las columnas es muy útil la
función janitor::clean_names() (en este caso ya la
aplicamos cuando generamos los datasets).
Esta librería también tiene otras funciones muy útiles, como
remove_empty() para eliminar filas o columnas vacías o
remove_constant() para eliminar filas o columnas
constantes.
Vamos a eliminar el sufijo del año en los nombres de las columnas porque sabemos que todos los datos son del mismo año:
df = df %>% rename_all(function(x) str_remove(x, "_2019"))
La librería skimr nos da un vistazo general de los
datos.
skim(df)
| Name | df |
| Number of rows | 5570 |
| Number of columns | 25 |
| _______________________ | |
| Column type frequency: | |
| character | 8 |
| numeric | 17 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| uf | 0 | 1 | 2 | 2 | 0 | 27 | 0 |
| municipio | 0 | 1 | 3 | 32 | 0 | 5570 | 0 |
| nome_regiao | 0 | 1 | 3 | 12 | 0 | 5 | 0 |
| nome_uf | 0 | 1 | 4 | 19 | 0 | 27 | 0 |
| semiarido | 0 | 1 | 3 | 3 | 0 | 2 | 0 |
| hierarquia | 0 | 1 | 9 | 18 | 0 | 5 | 0 |
| atividade_maior_vab | 0 | 1 | 10 | 84 | 0 | 10 | 0 |
| tipologia | 5 | 1 | 6 | 22 | 0 | 5 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| num_beneficiarios_mean | 0 | 1.00 | 2474.52 | 9208.35 | 2.00 | 346.71 | 947.50 | 2439.04 | 4.513838e+05 | ▇▁▁▁▁ |
| suma_valor | 0 | 1.00 | 5594117.72 | 18428002.57 | 3216.00 | 704285.50 | 2101369.00 | 5734926.75 | 8.539214e+08 | ▇▁▁▁▁ |
| codigo_uf | 0 | 1.00 | 32.38 | 9.83 | 11.00 | 25.00 | 31.00 | 41.00 | 5.300000e+01 | ▂▅▇▅▂ |
| codigo_municipio | 0 | 1.00 | 3253590.77 | 984910.34 | 1100015.00 | 2512125.75 | 3146280.00 | 4119189.50 | 5.300108e+06 | ▂▅▇▅▂ |
| vab_agro | 0 | 1.00 | 55783.48 | 101866.43 | 0.00 | 10113.36 | 25638.79 | 59948.02 | 1.575333e+06 | ▇▁▁▁▁ |
| vab_industria | 0 | 1.00 | 248797.85 | 1429430.02 | 449.30 | 4294.51 | 13302.82 | 72493.60 | 5.735987e+07 | ▇▁▁▁▁ |
| pib | 0 | 1.00 | 1326594.43 | 12696536.12 | 15515.27 | 87586.45 | 190970.75 | 510100.68 | 7.638060e+08 | ▇▁▁▁▁ |
| pib_percapita | 0 | 1.00 | 24546.81 | 25552.22 | 4482.85 | 10515.57 | 18182.28 | 29882.71 | 4.648835e+05 | ▇▁▁▁▁ |
| grado_urbanizacion | 5 | 1.00 | 48.28 | 34.28 | 0.00 | 0.00 | 54.78 | 77.68 | 9.967000e+01 | ▇▂▆▆▆ |
| area_geografica | 12 | 1.00 | 1526.55 | 5610.80 | 3.56 | 204.90 | 417.84 | 1025.39 | 1.595333e+05 | ▇▁▁▁▁ |
| despesa_ciencia_tecnologia | 289 | 0.95 | 63972.80 | 1805472.41 | 0.00 | 0.00 | 0.00 | 0.00 | 1.188419e+08 | ▇▁▁▁▁ |
| despesa_educacao_cultura | 289 | 0.95 | 33200752.21 | 214626480.31 | 0.00 | 5704397.12 | 11392275.52 | 24730333.47 | 1.358956e+10 | ▇▁▁▁▁ |
| despesa_pessoal_encargos_sociais | 286 | 0.95 | 64558552.22 | 463388560.08 | 278944.70 | 10610257.70 | 18681271.24 | 39789851.65 | 2.508161e+10 | ▇▁▁▁▁ |
| despesa_transporte | 289 | 0.95 | 2866400.07 | 72043623.38 | 0.00 | 5720.00 | 483191.64 | 1500556.47 | 5.171955e+09 | ▇▁▁▁▁ |
| exportacoes | 3379 | 0.39 | 102088499.29 | 484266387.22 | 1.00 | 350150.50 | 4824744.00 | 41587247.00 | 1.116350e+10 | ▇▁▁▁▁ |
| importacoes | 3553 | 0.36 | 87857804.51 | 514210623.13 | 48.00 | 84013.00 | 1237977.00 | 12267004.00 | 1.012875e+10 | ▇▁▁▁▁ |
| populacao | 13 | 1.00 | 37600.16 | 221267.52 | 781.00 | 5446.00 | 11638.00 | 25501.00 | 1.225202e+07 | ▇▁▁▁▁ |
skim_eda = partition(skim(df))
names(skim_eda)
## [1] "character" "numeric"
Pueden revisar la
vignette de skimr para conocer más opciones.
Calculamos el % de datos faltantes por variable:
map_dbl(df, function(x) mean(is.na(x)) * 100) %>%
sort(decreasing=T) %>%
as.data.frame()
## .
## importacoes 63.78815081
## exportacoes 60.66427289
## despesa_ciencia_tecnologia 5.18850987
## despesa_educacao_cultura 5.18850987
## despesa_transporte 5.18850987
## despesa_pessoal_encargos_sociais 5.13464991
## populacao 0.23339318
## area_geografica 0.21543986
## grado_urbanizacion 0.08976661
## tipologia 0.08976661
## uf 0.00000000
## num_beneficiarios_mean 0.00000000
## suma_valor 0.00000000
## municipio 0.00000000
## nome_regiao 0.00000000
## codigo_uf 0.00000000
## nome_uf 0.00000000
## codigo_municipio 0.00000000
## semiarido 0.00000000
## hierarquia 0.00000000
## vab_agro 0.00000000
## vab_industria 0.00000000
## pib 0.00000000
## pib_percapita 0.00000000
## atividade_maior_vab 0.00000000
y por municipio:
na_by_row = df %>% apply(1, function(x) mean(is.na(x)) * 100)
names(na_by_row) = df$municipio
na_by_row %>% sort(decreasing=T) %>% head(10)
## balneario_rincao mojui_dos_campos paraiso_das_aguas pescaria_brava
## 40 40 40 40
## pinto_bandeira campo_grande_3 conchas entre_rios_2
## 40 32 32 32
## muribeca nazaria
## 32 32
mean(na_by_row > 0)
## [1] 0.7204668
También podemos recurrir a skimr:
skim_eda$numeric %>%
select(skim_variable, n_missing) %>%
filter(n_missing > 0)
Variable type: numeric
| skim_variable | n_missing |
|---|---|
| grado_urbanizacion | 5 |
| area_geografica | 12 |
| despesa_ciencia_tecnologia | 289 |
| despesa_educacao_cultura | 289 |
| despesa_pessoal_encargos_sociais | 286 |
| despesa_transporte | 289 |
| exportacoes | 3379 |
| importacoes | 3553 |
| populacao | 13 |
skim_eda$character %>%
select(skim_variable, n_missing) %>%
filter(n_missing > 0)
Variable type: character
| skim_variable | n_missing |
|---|---|
| tipologia | 5 |
En un caso real hay que ser muy cuidadoso con los datos faltantes. El tratamiento que les demos depende del tipo de modelo predictivo que usemos: en algunos casos vamos a necesitar imputarlos mientras que ciertas metodologías admiten NAs.
Durante el EDA tomaremos las decisiones que más nos sirvan para entender los datos y comunicar los resultados. Para esto es fundamental preguntarse por la interpretación de los missing variable por variable.
En este ejemplo vamos a suponer que los NA en las variables monetarias indican que el monto es igual a 0.
vars_monetarias = c(
"despesa_ciencia_tecnologia", "despesa_educacao_cultura",
"despesa_pessoal_encargos_sociais", "despesa_transporte",
"exportacoes", "importacoes", "pib", "pib_percapita",
"suma_valor", "vab_agro", "vab_industria"
)
df = df %>%
mutate_at(vars_monetarias, function(x) ifelse(is.na(x), 0, x))
Además vemos que los datos de población con NA se pueden reconstruir con el PIB y el PIB per cápita.
df = df %>%
mutate(
populacao = ifelse(
is.na(populacao), pib / pib_percapita * 1000, populacao)
,populacao = round(populacao, 0)
)
Para explorar los datos puede ser útil crear nuevas columnas o transformar las que tenemos.
df = df %>%
mutate(
densidad = populacao / area_geografica
,suma_perbenef = suma_valor / num_beneficiarios_mean
)
df = df %>%
mutate_at(
vars(
vab_agro, vab_industria, despesa_ciencia_tecnologia
,despesa_pessoal_encargos_sociais, despesa_educacao_cultura
,despesa_transporte, exportacoes, importacoes
,num_beneficiarios_mean, suma_valor
), list(pc=function(x) x / .$populacao)
)
# recodificar variables
df = df %>%
mutate_at(vars(codigo_municipio, codigo_uf), as.character)
# categorías nuevas en función de valores de otras variables
df = df %>%
mutate(
semiarido = ifelse(semiarido == "Sim", T, F),
regiao = case_when(
nome_regiao == "Norte" ~ "N"
,nome_regiao == "Nordeste" ~ "NE"
,nome_regiao == "Centro-oeste" ~ "CO"
,nome_regiao == "Sudeste" ~ "SE"
,nome_regiao == "Sul" ~ "S"
,TRUE ~ "error"
)
)
Podemos usar forcats::fct_lump() para agrupar categorías
poco frecuentes.
df$atividade_maior_vab %>% unique()
## [1] "Demais serviços"
## [2] "Administração, defesa, educação e saúde públicas e seguridade social"
## [3] "Agricultura, inclusive apoio à agricultura e a pós colheita"
## [4] "Indústrias de transformação"
## [5] "Pecuária, inclusive apoio à pecuária"
## [6] "Eletricidade e gás, água, esgoto, atividades de gestão de resíduos e descontaminação"
## [7] "Comércio e reparação de veículos automotores e motocicletas"
## [8] "Produção florestal, pesca e aquicultura"
## [9] "Indústrias extrativas"
## [10] "Construção"
df %>%
mutate(
atividade = forcats::fct_lump(
atividade_maior_vab, prop=0.1, other_level="Otros")
) %>%
pull(atividade) %>%
unique()
## [1] Demais serviços
## [2] Administração, defesa, educação e saúde públicas e seguridade social
## [3] Agricultura, inclusive apoio à agricultura e a pós colheita
## [4] Otros
## 4 Levels: Administração, defesa, educação e saúde públicas e seguridade social ...
A veces es útil separar el dataset según el tipo de variables una vez que terminamos el preprocesamiento.
df_num = df %>% select_if(is.numeric)
df_cat = df %>% select_if(function(x) !is.numeric(x))
Usamos histogramas y densidades para visualizar la distribución.
ggplot(df, aes(x=pib_percapita)) +
geom_histogram(alpha=0.5) +
NULL
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(df, aes(x=log10(pib_percapita))) +
geom_histogram(alpha=0.5) +
NULL
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La escala logarítmica es una escala no lineal en la que moverse una unidad implica multiplicarse por las base del logaritmo (2, 10 o e, por ejemplo). \(\text{log}_{10}\) mide, por ejemplo, órdenes de magnitud (\(x\) difiere en un orden de magnitud de \(y\) si \(x\) es aproximadamente 10 veces mayor que \(y\)).
Entonces, la transformación logarítmica es útil cuando la variabilidad relevante está en el orden de magnitud, no en la escala original. Esto es común en el caso de variables monetarias, cuando las variaciones relevantes son porcentuales y no nominales.
Otra manera de visualizar lo mismo:
ggplot(df, aes(x=pib_percapita)) +
geom_histogram(alpha=0.5) +
scale_x_log10(
breaks = trans_breaks("log10", function(x) 10^x),
labels=trans_format("log10", math_format(10^.x))
) +
NULL
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Podemos comparar la distribución de una variable entre grupos.
ggplot(df, aes(x=log(pib_percapita), fill=semiarido)) +
geom_histogram(alpha=0.5) +
NULL
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# cuando usamos density normalizamos la distribucion en cada grupo
ggplot(df, aes(x=log(pib_percapita), fill=semiarido)) +
geom_density(alpha=0.5, adjust=2) +
NULL
# adjust determina el nivel de suavidad
# adjust=2 --> bandas de 2 veces el ancho default
Con la librería ggridges podemos hacer algunos gráficos
lindos:
ggplot(df, aes(x=log(pib_percapita), y=tipologia, fill =..x..)) +
ggridges::geom_density_ridges_gradient(scale=1) +
scale_fill_viridis_c() +
NULL
## Picking joint bandwidth of 0.188
ggplot(df, aes(x=log(pib_percapita), y=tipologia, fill=tipologia)) +
ggridges::stat_density_ridges(
quantile_lines=T, quantiles=c(.25,.5,.75), alpha=0.7) +
NULL
## Picking joint bandwidth of 0.188
# a veces es útil visualizar la cantidad de obs al mismo tiempo!
ggplot(df, aes(x=log(pib_percapita), y=tipologia, fill=tipologia)) +
ggridges::geom_density_ridges(jittered_points=T, point_size=0.5) +
NULL
## Picking joint bandwidth of 0.188
También podemos usar boxplots:
Fuente: https://r4ds.had.co.nz/
ggplot(df, aes(
x=fct_reorder(nome_regiao, pib_percapita, .fun="median"),
y=pib_percapita, fill=nome_regiao)) +
geom_boxplot() +
scale_y_continuous(trans='log10') +
# facet_wrap(~ tipologia) +
theme_minimal() +
NULL
# eliminando los outliers, invirtiendo las coordenadas, y mostrando los puntos:
ggplot(df, aes(
x=fct_reorder(hierarquia, pib_percapita, .fun="median"),
y=pib_percapita, fill=hierarquia)
) +
geom_jitter(aes(color=hierarquia), alpha=0.2, size=0.5) +
gg.layers::geom_boxplot2(alpha=0.5) +
scale_y_continuous(trans='log10') +
# facet_wrap(vars(nome_regiao)) +
coord_flip() +
theme_minimal() +
theme(legend.position="none") +
NULL
Podemos calcular estadísticos globales o por grupo:
quantile(df_num$pib_percapita, seq(0,1,0.1))
## 0% 10% 20% 30% 40% 50% 60%
## 4482.850 8124.319 9620.964 11579.717 14544.326 18182.285 22222.446
## 70% 80% 90% 100%
## 26959.432 33584.336 45617.949 464883.490
# para todas las variables:
list_cuantiles = df_num %>%
map(function(x) quantile(x, seq(0, 1, 0.25), na.rm=T) )
list_cuantiles %>%
as.data.frame() %>%
select(num_beneficiarios_mean, suma_valor, pib)
## num_beneficiarios_mean suma_valor pib
## 0% 2.0000 3216.0 15515.27
## 25% 346.7083 704285.5 87586.45
## 50% 947.5000 2101369.0 190970.75
## 75% 2439.0417 5734926.8 510100.68
## 100% 451383.8333 853921396.0 763805984.80
¡También podemos revisar los resultados de skimr!
# por grupo:
df %>%
group_by(nome_regiao) %>%
summarise(
n = n(),
expo_pc = mean(exportacoes_pc, na.rm=T)
,expo_pc2 = weighted.mean(
exportacoes_pc, w=populacao, na.rm=T)
,popM = sum(populacao, na.rm=T) / 1e6
,area = sum(area_geografica, na.rm=T)
)
## # A tibble: 5 x 6
## nome_regiao n expo_pc expo_pc2 popM area
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Centro-oeste 467 2132. 1601. 16.1 1557857.
## 2 Nordeste 1794 159. 291. 60.6 1630331.
## 3 Norte 450 704. 1214. 16.7 3796198.
## 4 Sudeste 1668 948. 1302. 87.2 919671.
## 5 Sul 1191 669. 1526. 30.5 580495.
# Cuales son las expo. per capita de cada region?
En general usamos diagramas de dispersión (scatter plots):
ggplot(df, aes(x=suma_perbenef, y=pib_percapita)) +
geom_point() +
NULL
ggplot(df, aes(x=log10(suma_perbenef), y=log10(pib_percapita))) +
geom_point(alpha=0.3, size=0.5) +
geom_smooth() + # smoother / suavizado
NULL
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(df,
aes(x=log(suma_valor_pc), y=log(pib_percapita),
color=semiarido)) +
geom_point(size=0.5, alpha=0.5) +
geom_smooth() + # smoother / suavizado
NULL
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Una relación lineal en un gráfico log-log implica elasticidad constante: la pendiente de un ajuste lineal indica aproximadamente cuánto crece en promedio porcentualmente una variable cuando la otra crece 1%.
# a veces es mejor un plot por grupo:
ggplot(df, aes(x=num_beneficiarios_mean_pc, y=log10(pib_percapita))) +
geom_point(size=.7, alpha=0.3) +
geom_smooth(se=F, method="loess", size=.8) +
facet_wrap(vars(nome_regiao), scales="free_y") +
theme_light() +
NULL
## `geom_smooth()` using formula 'y ~ x'
Cuando N es grande, podemos usar alguna de las siguientes opciones:
# partir en cuadrados:
ggplot(df) +
geom_bin2d(aes(log10(suma_valor_pc), log10(pib_percapita))) +
scale_fill_viridis_c()
# partir en hexagonos:
# install.packages("hexbin")
ggplot(df) +
geom_hex(aes(log10(suma_valor_pc), log10(pib_percapita))) +
scale_fill_viridis_c()
# usar boxplots:
ggplot(df, aes(log10(suma_valor_pc), log10(pib_percapita))) +
geom_boxplot(aes(group=cut_number(log10(suma_valor_pc), 10)))
# todos los grupos tienen igual N
# despegar puntitos pegados (con cuidado, p/variables discretas):
# ggplot(df) +
# geom_jitter(aes(log10(suma_valor_pc), log10(pib_percapita)))
Con un correlograma visualizamos la correlación entre todas las variables cuantitativas:
GGally::ggcorr(
df_num, method=c("pairwise","spearman"),
label=T, hjust=1, label_size=2, layout.exp=10, size=3)
Las correlaciones de a pares pueden ser útiles para detectar variables redundantes:
cor_matrix = cor(df_num, method="spearman", use="pairwise")
cor_matrix[upper.tri(cor_matrix, diag=T)] = NA
df_cor = cor_matrix %>% as.table() %>% as.data.frame()
df_cor %>%
rename(corr = Freq) %>%
filter(!is.na(corr) & Var1 != Var2) %>%
arrange(-abs(corr)) %>%
head(10) %>%
knitr::kable() %>%
kableExtra::kable_styling()
| Var1 | Var2 | corr |
|---|---|---|
| despesa_ciencia_tecnologia_pc | despesa_ciencia_tecnologia | 0.9998316 |
| importacoes_pc | importacoes | 0.9961378 |
| exportacoes_pc | exportacoes | 0.9932097 |
| suma_valor | num_beneficiarios_mean | 0.9908110 |
| suma_valor_pc | num_beneficiarios_mean_pc | 0.9829813 |
| despesa_pessoal_encargos_sociais | despesa_educacao_cultura | 0.9586380 |
| pib | vab_industria | 0.9265616 |
| despesa_transporte_pc | despesa_transporte | 0.8852189 |
| populacao | despesa_educacao_cultura | 0.8561423 |
| populacao | despesa_pessoal_encargos_sociais | 0.8551407 |
Métodos de manejo de NA en
cor():
"everything": usa todas las observaciones (si hay NA en un par de variables, arroja NA)"pairwise": usa las observaciones completas de a pares (descarta las observaciones con algún NA de a pares)"complete.obs": usa las observaciones completas del data.frame (descarta las observaciones con algún NA considerando todas las filas)
En la clase que viene vamos a ver cómo medir asociaciones entre distintos tipos de variables, y cómo esto nos puede ayudar para guiar el EDA.
Revisamos la distribución de frecuencias de las variables categóricas.
table(df_cat$tipologia)
##
## IntermediarioAdjacente IntermediarioRemoto RuralAdjacente
## 686 60 3040
## RuralRemoto Urbano
## 323 1456
table(df_cat$tipologia, df_cat$semiarido)
##
## FALSE TRUE
## IntermediarioAdjacente 525 161
## IntermediarioRemoto 54 6
## RuralAdjacente 2162 878
## RuralRemoto 240 83
## Urbano 1322 134
# para todas las variables:
tabs = df_cat %>% map(function(x) table(x, useNA="ifany"))
tabs$hierarquia
## x
## Capital Regional Centro de Zona Centro Local Centro Subregional
## 189 561 4479 164
## Metrópole
## 177
tabs$tipologia
## x
## IntermediarioAdjacente IntermediarioRemoto RuralAdjacente
## 686 60 3040
## RuralRemoto Urbano <NA>
## 323 1456 5
Podemos hacer tablas más informativas con
janitor::tabyl().
df_cat %>%
tabyl(tipologia, semiarido) %>%
adorn_totals() %>%
adorn_percentages("row") %>%
adorn_pct_formatting() %>%
adorn_ns()
## tipologia FALSE TRUE
## IntermediarioAdjacente 76.5% (525) 23.5% (161)
## IntermediarioRemoto 90.0% (54) 10.0% (6)
## RuralAdjacente 71.1% (2162) 28.9% (878)
## RuralRemoto 74.3% (240) 25.7% (83)
## Urbano 90.8% (1322) 9.2% (134)
## <NA> 100.0% (5) 0.0% (0)
## Total 77.3% (4308) 22.7% (1262)
Cortamos el target numérico en bins y hacemos tablas de frecuencias:
tablas_vs_target = function(target, df_cat, variable, g=5) {
df_cat$target_ = Hmisc::cut2(target, g=g)
tab = df_cat %>%
tabyl(target_, !!sym(variable)) %>%
adorn_totals() %>%
adorn_percentages("col") %>%
adorn_pct_formatting() %>%
adorn_ns()
return(tab)
}
tablas = list()
for (col in names(df_cat)) {
tablas[[col]] = tablas_vs_target(
df$pib_percapita, df_cat, col)
}
tablas$hierarquia
## target_ Capital Regional Centro de Zona Centro Local
## [ 4483, 9622) 2.1% (4) 9.8% (55) 23.4% (1046)
## [ 9622, 14545) 5.8% (11) 19.4% (109) 21.3% (953)
## [14545, 22225) 19.0% (36) 16.4% (92) 20.0% (894)
## [22225, 33591) 25.9% (49) 27.8% (156) 18.8% (841)
## [33591,464883] 47.1% (89) 26.6% (149) 16.6% (745)
## Total 100.0% (189) 100.0% (561) 100.0% (4479)
## Centro Subregional Metrópole
## 3.0% (5) 2.3% (4)
## 11.6% (19) 12.4% (22)
## 24.4% (40) 29.4% (52)
## 26.2% (43) 14.1% (25)
## 34.8% (57) 41.8% (74)
## 100.0% (164) 100.0% (177)
tablas$semiarido
## target_ FALSE TRUE
## [ 4483, 9622) 9.9% (425) 54.6% (689)
## [ 9622, 14545) 16.1% (695) 33.2% (419)
## [14545, 22225) 23.4% (1009) 8.3% (105)
## [22225, 33591) 25.2% (1085) 2.3% (29)
## [33591,464883] 25.4% (1094) 1.6% (20)
## Total 100.0% (4308) 100.0% (1262)
Podemos hacer tablas más lindas con knitr::kable() y
kableExtra::kable_styling().
tablas$hierarquia %>%
knitr::kable() %>%
kableExtra::kable_styling(
bootstrap_options=c("striped", "hover"), font_size=12
)
| target_ | Capital Regional | Centro de Zona | Centro Local | Centro Subregional | Metrópole |
|---|---|---|---|---|---|
| [ 4483, 9622) | 2.1% (4) | 9.8% (55) | 23.4% (1046) | 3.0% (5) | 2.3% (4) |
| [ 9622, 14545) | 5.8% (11) | 19.4% (109) | 21.3% (953) | 11.6% (19) | 12.4% (22) |
| [14545, 22225) | 19.0% (36) | 16.4% (92) | 20.0% (894) | 24.4% (40) | 29.4% (52) |
| [22225, 33591) | 25.9% (49) | 27.8% (156) | 18.8% (841) | 26.2% (43) | 14.1% (25) |
| [33591,464883] | 47.1% (89) | 26.6% (149) | 16.6% (745) | 34.8% (57) | 41.8% (74) |
| Total | 100.0% (189) | 100.0% (561) | 100.0% (4479) | 100.0% (164) | 100.0% (177) |
Y aún más lindas con gt
o gtsummary
:)
En general usamos gráficos de barras:
ggplot(df) +
geom_bar(aes(x=nome_regiao)) +
NULL
Vamos a mejorarlo un poco:
# ordenando por N
ggplot(df) +
geom_bar(aes(x=fct_infreq(nome_regiao))) +
NULL
Agregando labels:
# La version mas general: DF agrupado + geom_col()
# (en lugar de geom_bar())
df_tmp = df %>%
group_by(nome_regiao) %>%
summarise(
n = n()
) %>%
mutate(
porcentaje = n / sum(n) * 100,
nome_regiao = fct_reorder(nome_regiao, n)
)
ggplot(df_tmp) +
geom_col(aes(x=nome_regiao, y=porcentaje), alpha=0.5) +
geom_text(aes(
x=nome_regiao, y=porcentaje, label=glue("{n}\n({round(porcentaje,2)}%)")),
color="black", size=3) +
theme_light() +
labs(x="Region", y="N") +
NULL
Nota: a veces con una buena tabla alcanza!
Podemos analizar la relación con alguna variable numérica usando alguna medida resumen como la media:
ggplot(df) +
geom_bar(aes(x=nome_regiao, y=pib_percapita),
fun="mean", stat="summary", alpha=.5) +
theme_minimal() +
NULL
Los puntos suelen usarse para representar valores donde no es necesario incluir el 0.
ggplot(df, aes(x=nome_regiao, y=pib_percapita)) +
geom_point(fun="mean", stat="summary") +
# stat_summary(fun = "mean", geom = "point") +
geom_line(aes(group=1), fun="mean", stat="summary") +
theme_minimal() +
NULL
Ordenando los niveles según los valores de otra variable:
# resumir los datos y luego usar geom_col():
df_tmp = df %>%
group_by(nome_uf) %>%
summarise_at(vars(populacao), sum)
ggplot(df_tmp) +
geom_col(
aes(x=populacao, y=fct_reorder(as.factor(nome_uf), populacao)),
alpha=0.5) +
labs(y=NULL) +
scale_x_continuous(labels=scales::comma) +
theme_light() +
NULL
# # La otra manera:
# ggplot(df) +
# geom_bar(
# aes(x=populacao, y=fct_reorder(as.factor(nome_uf), populacao, .fun=sum)),
# fun="sum", stat="summary", alpha=0.5) +
# labs(y=NULL) +
# scale_x_continuous(labels = scales::comma) +
# theme_classic() +
# NULL
También podemos hacer un boxplot por grupo, como vimos antes.
Vamos a suponer momentáneamente que nuestro target es binario: PIB bajo (0) o PIB alto (1). ¿Qué tablas y visualizaciones nos pueden servir como EDA en este escenario?
df = df %>%
mutate(target_pib = ifelse(pib_percapita > mean(pib_percapita), 1, 0))
# TODO
Con plotly::ggplotly() podemos hacer interactivo
cualquier gráfico de ggplot2.
g = ggplot(
df
,aes(x=log(pib_percapita+1), y=log(exportacoes+1)
,color=semiarido, label=municipio, label2=nome_regiao)) +
geom_point(alpha=0.8, size=0.8) +
theme_minimal() +
NULL
plotly::ggplotly(g)
Con GGally::ggpairs() podemos analizar relaciones entre
variables de múltiples tipos.
df_reduced = df %>%
select(
nome_regiao, suma_perbenef, densidad,
despesa_ciencia_tecnologia_pc, pib_percapita, semiarido
) %>%
mutate_if(is.numeric, function(x) log10(x+1))
GGally::ggpairs(
df_reduced,
columns = c("suma_perbenef", "densidad",
"despesa_ciencia_tecnologia_pc", "pib_percapita"),
mapping = aes(color=semiarido),
lower = list(continuous = GGally::wrap("points", alpha=0.3, size=0.1))
)
¿Cómo generamos todos los gráficos bivariados con respecto a la variable target?
# generamos log de todas las vars
df_num_expanded = df_num %>%
mutate_all(list("log" = function(x) ifelse(x==0, 0, sign(x)*log10(abs(x))) ))
# Medimos corr en logs vs el target (equivale a spearman!)
corr_mat = cor(
df_num_expanded %>% select(ends_with("_log")) %>% select(-pib_percapita_log),
df_num_expanded %>% select(pib_percapita_log),
method="pearson", use="pairwise")
df_corr = corr_mat %>% as.data.frame()
names(df_corr) = "correlation"
df_corr = df_corr %>%
arrange(-abs(correlation)) %>%
rownames_to_column("variable")
# plot function
plt_corr = function(df, variable) {
p = ggplot(df, aes(x=!!sym(variable), y=pib_percapita_log)) +
geom_jitter(size=.3, alpha=.2) +
geom_smooth(se=T, method="loess", size=.8, formula=y~x) +
labs(title=variable) +
theme_minimal() +
NULL
return(p)
}
# lista de plots
variables = df_corr$variable
list_plots = map(variables, function(x) plt_corr(df_num_expanded, x))
# save plots
ggsave("output/bivariado_pib.pdf",width=6, height=6,
gridExtra::marrangeGrob(grobs=list_plots, nrow=2, ncol=1))
Descargamos los shapefiles de IBGE.
Shapefile es un formato de archivo con información geográfica. Indican el sistemas de coordenadas de referencia y proyección cartográfica usados. Contiene puntos (posiciones), líneas (recorridos) y polígonos (superficies) de alguna parte del mundo.
Los podemos leer con sf::read_sf() y hacer un
join con nuestros datos.
sf_mun = sf::read_sf(
"data/raw/ibge/br_municipios_20200807/BR_Municipios_2019.shp")
sf_uf = sf::read_sf(
"data/raw/ibge/br_unidades_da_federacao/BR_UF_2019.shp")
names(sf_mun)
## [1] "CD_MUN" "NM_MUN" "SIGLA_UF" "AREA_KM2" "geometry"
dim(sf_mun)
## [1] 5572 5
sf_mun = sf_mun %>%
left_join(df, by=c("CD_MUN"="codigo_municipio"))
# municipios faltantes en nuestros datos
sf_mun %>%
filter(is.na(municipio)) %>%
select(CD_MUN, NM_MUN, SIGLA_UF, AREA_KM2)
## Simple feature collection with 2 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -53.53194 ymin: -33.61715 xmax: -50.53778 ymax: -30.00031
## Geodetic CRS: SIRGAS 2000
## # A tibble: 2 x 5
## CD_MUN NM_MUN SIGLA_UF AREA_KM2 geometry
## <chr> <chr> <chr> <dbl> <MULTIPOLYGON [°]>
## 1 4300001 Lagoa Mirim RS 2872. (((-53.52638 -33.61715, -53.5284 -3~
## 2 4300002 Lagoa dos Patos RS 10197. (((-52.09772 -32.16174, -52.09888 -~
Creamos mapas coropléticos con geom_sf(). En estos mapas
las áreas se sombrean de distintos colores según alguna
característica.
sf_mun_tmp = sf_mun %>% filter(nome_uf == "Sergipe")
sf_uf_tmp = sf_uf %>% filter(NM_UF == "Sergipe")
ggplot() +
geom_sf(data=sf_mun_tmp, aes(fill=pib_percapita), color="black", size=0.5) +
geom_sf(data=sf_uf_tmp, fill=NA, color="black", size=0.8) +
scale_fill_viridis_c(trans="log10", na.value="gray60") +
# coord_sf(xlim=c(-75,-55), ylim=c(-10, 3)) +
labs(fill="PIB per capita") +
NULL
g = ggplot() +
geom_sf(data=sf_mun, aes(fill=pib_percapita), color=NA) +
geom_sf(data=sf_uf, fill=NA, color="black", size=0.5) +
scale_fill_viridis_c(trans="log10", na.value="gray60") +
theme_minimal() +
guides(scale ="none") +
theme(legend.position=c(0.15,0.2)) +
labs(fill="PIB per capita") +
NULL
ggsave("output/mapa_pibpercapita.png", g, width=10, height=8)
Si tuviéramos lat y long en nuestros datos podríamos inferir el mapa de esta manera:
df_geo = df %>%
st_as_sf(coords=c("longitud", "latitud"), crs=4326)
ggplot() + geom_sf(df_geo)